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PURPOSE:  Level  set  methods  are  often  used  to  capture  interface  behavior  in  two-phase,  in¬ 
compressible  flow  models.  While  level  set  techniques  for  structured  computational  grids  have  been 
widely  investigated,  approaches  for  unstructured  meshes  are  less  mature.  This  report  details  the 
formulation  and  implementation  of  a  discontinuous  Galerkin-based  approach  that  is  suitable  for 
unstructured  meshes  and  offers  potential  gains  in  accuracy  and  efficiency  over  more  traditional 
level  set  techniques. 

INTRODUCTION:  Flow  of  air  and  water  around  solid  objects  can  be  modeled  by  assuming 
the  phases  are  separated  by  sharp  interfaces  and  that  each  fluid  subdomain  is  governed  by  the 
Navier-Stokes  equations.  The  resulting  model  requires  resolution  of  the  flow  field  and  the  location 
and  evolution  of  the  interfaces.  There  are  many  techniques  available  for  approximating  flow  in 
each  subdomain  and  many  ways  to  resolve  the  interfaces.  Level  set  methods  represent  one  such 
class  of  techniques  for  capturing  interface  behavior  that  has  been  applied  successfully  to  problems 
in  fields  from  computational  geometry  to  fluid  mechanics  (Osher  and  Fedkiw  2001,  Sethian  2001). 
Among  other  things,  the  appeal  of  level  set  methods  can  be  attributed  to  the  generality  of  their 
formulation,  ability  to  resolve  interfaces  accurately  while  allowing  for  topological  changes,  and 
relatively  straightforward  extension  to  higher  dimensional  problems  (Sethian  1999). 

Level  set  methods  for  structured  grids  are  fairly  mature  (Sethian  2001,  Osher  and  Fedkiw  2001). 
Although  extensions  for  general  interface  propagation  problems  as  well  as  two-phase  flow  have 
been  considered  (Barth  and  Sethian  1998,  Sethian  and  Vladimirsky  2000,  Nagrath  et  al.  2005, 
Smolianski  2005),  approaches  for  unstructured  meshes  are  less  mature.  With  this  in  mind,  we  are 
interested  in  the  design  and  implementation  of  level  set  techniques  for  air/water  flow  that  naturally 
apply  to  unstructured  tetrahedral  meshes  and  are  well  suited  to  distributed  computing  architectures. 

Standard  level  set  methodology  builds  upon  an  accurate  discretization  for  a  linear  hyperbolic  partial 
differential  equation  (PDE)  or  an  equivalent  Hamilton- Jacobi  formulation  (Sethian  2001).  It  also 
relies  heavily  on  an  effective  reinitialization  technique  that  must  be  employed  at  intermediate  times 
to  recover  the  accuracy  of  the  level  set  representation  of  the  fluid  front  (Sussman  and  Fatemi  1999). 

In  this  report  we  focus  on  the  implementation  of  an  alternative  strategy  that  employs  a  Runge-Kutta 
discontinuous  Galerkin  (RKDG)  discretization  (Cockbum  and  Shu  2001)  and  exploits  the  fluid 
incompressibility  assumption  (Marchandise  et  al.  2006).  With  sufficiently  high-order  approxima¬ 
tions,  this  method  potentially  offers  accurate  solutions  with  good  mass  conservation,  while  obviat¬ 
ing  the  reinitialization  step.  It  also  lends  itself  to  fully  explicit  time  integration  and  a  quadrature- 
free  implementation  that  is  readily  parallelizable  (Atkins  and  Shu  1996,  Baggag  et  al.  1999). 

FORMULATION:  We  begin  with  a  physical  domain  Q  in  which  an  interface  F(/:)  evolves  over 
a  time  interval  [0,  T\.  To  characterize  I  ff),  we  adopt  a  level  set  formulation  and  define  implicitly 
a  function  </>(x,  t )  such  that  0(x,  t)  =  0  corresponds  to  T(t).  In  this  case,  the  propagation  of  T(t) 
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with  normal  speed  un  can  be  expressed  as  (Osher  and  Fedkiw  2001,  Sethian  2001) 


(1) 


with  the  initial  data,  0(x,  0)  =  0°(x),  chosen  so  that  0°(x)  =  0  on  T(0).  We  further  require  that 
(j>°(x)  be  a  signed  distance  (i.e.,  0°(x)  =  ±d  where  d  is  the  distance  to  T),  although  this  is  not 
strictly  necessary  (Olsson  and  Kreiss  2005). 

Successful  level  set  approximations  require  accurate  solution  of  Equation  1  and  a  suitable  velocity, 
u,  defined  throughout  Q  that  gives  the  correct  front  propagation  speed,  un,  on  T  (Sethian  2001).  In 
addition,  these  methods  typically  require  an  efficient  approach  for  initializing  0(x,  t  =  trn )  at  given 
time  instances,  tm,  that  ensures  0(x,  tm )  is  sufficiently  smooth  over  f)  but  yet  maintains  0(x,  tm )  = 
0  on  T(tm)  (Sussman  and  Fatemi  1999).  As  mentioned  in  the  introduction,  the  emphasis  here  is  on 
one  configuration  of  discrete  approximations  and  solution  algorithms  that  attempts  to  address  the 
above  requirements  and  may  hold  promise  for  simulating  multiphase,  incompressible  fluid  flow. 

Conservative  level  set  equation  formulation:  The  approaches  considered  rely 

on  discontinuous  Galerkin  (DG)  spatial  discretizations  to  obtain  accurate  numerical  solutions  for 
0  on  unstructured  meshes.  DG  approximations  for  0  can  be  obtained  in  at  least  two  ways.  Equa¬ 
tion  1  can  be  solved  directly  using  a  general  RKDG  formulation  for  Hamilton-Jacobi  equations 
(Hu  and  Shu  1999,  Li  and  Shu  2005).  Or,  we  can  take  advantage  of  the  fact  that  our  area  of  interest 
is  incompressible,  multiphase  flow  in  order  to  apply  methods  for  solving  conservative,  linear  ad- 
vection  problems.  We  follow  Marchandise  et  al.  (2006)  and  adopt  the  latter  approach  here,  since  it 
allows  us  to  use  tools  we  have  previously  developed  for  solving  PDEs  with  DG  methods  (Li  et  al. 
2007). 

First,  Equation  1  is  equivalent  to  solving 


(2) 


for  an  appropriately  defined  u.  In  the  case  of  two-phase  flow,  we  can  identify  u  with  the  fluid 
velocity  so  that  u  ■  n  =  un,  where  n  =  V0/||  V0||  is  the  unit  normal  to  the  interface  T  (Chang 
et  al.  1996).  Equation  2  can  then  be  rearranged  as 


—  +  V  ■  (U0)  =  0V  ■  u 


(3) 


which  gives 


(4) 


when  the  fluid  is  incompressible  (i.e.,  V  •  u  =  0)  (Marchandise  et  al.  2006). 

Element  weak  formulation:  Equation  4  is  just  the  linear  advection  equation  in  conser¬ 
vative  form,  and  a  weak  formulation  follows  in  a  standard  way  (Cockbum  and  Shu  2001).  Given  a 
triangulation,  A4h,  of  Q  and  an  element  8  €  A4h,  an  approximate  solution  is  sought  in  the  space 


vh  =  {vh  e  L°°(D)  :  vh\e  G  Vh(£),V£  G  Mh } 


(5) 
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where  we  denote  the  local  discrete  test  and  trial  space  as  Vh(S).  Multiplying  Equation  4  by  a  test 
function  and  integrating  by  parts  over  an  element  £,  we  have 


vh-^-dx=  /  (f)u-Vvhdx-  /  vh(f)u-n£ds,  Vvh  eVh(£) 


dt 


'as 


where  n£  is  the  unit  outer  normal  on  £.  Inserting  a  trial  solution  0^  e  Vh{£)  gives 
[  vh^-dx=  [  (j)hu  ■  Vvh  dx  -  [  vhcj)hu- n£ds,  Vvh  eVh(£) 


dt 


i as 


(6) 


(7) 


Since  the  underlying  spaces  are  discontinuous,  the  flux  on  internal  element  boundaries  is  multiply 
defined.  We  then  replace  the  outer  flux  in  the  last  term  of  the  right-hand  side  of  Equation  7  with  a 
numerical,  upwinded  flux 


where 


vh- 


94>h 

dt 


dx=  0/iU  •  V'i'h  dx  — 


'as 


vh(j)lpu  ■  n£  ds,  Vvh  e  Vh(£) 


cj)up 

<t>~ 

0+ 


J  0  u  •  n£  >  0 
[  0+  otherwise 

lim  0(x  +  en^,  t ) 

e— >0“ 

lim  0(x  +  en^,f) 

e— >0+ 


(8) 


(9) 


That  is,  oup  is  simply  the  value  of  0  taken  from  the  upwind  element  at  an  element  interface. 

DISCRETE  APPROXIMATION:  Equation  8  represents  a  semi-discrete  system  on  each 
element  of  A4h  with  coupling  across  elements  introduced  through  element  boundary  fluxes.  In  the 
following,  we  detail  an  RKDG  discrete  approximation  for  Equation  8  on  affine,  simplicial  meshes, 
which  is  appealing  in  its  simplicity. 

Finite  element  approximation:  We  assume  that  there  is  a  consistently  defined,  affine 
mapping  F  from  a  reference  element,  £,  for  all  £  6  A4h,  where  £  is  the  reference  simplex  in  Wl, 
d  =  1,2,3  (see  Figure  1).  A  local  basis  for  £  is  denoted  {A^}  for  i  —  1, . . . ,  np,  where 

N,  =  Ni  o  F"1 

and  {Ni}  is  a  basis  for  Pk(£),  the  space  of  polynomials  on  £  of  degree  k  or  lower.  nv  is  the 
dimension  of  Pk{£).  To  define  the  local  shape  functions,  {Ni},  we  consider  the  standard  nodal 
Lagrangian  polynomials  on  £  (see  Figure  2  for  the  Pi{£)  case).  We  label  the  corresponding  set  of 
nodal  locations  {pt}  with  iVj(pj)  =  dt],  for  i,  j  =  1, . . . ,  np.  In  general,  we  follow  the  convention 
of  using  a  hat  Q  to  signify  quantities  associated  with  a  reference  element  and  write  a  trial  solution 
as  4>h  =  1 0J7Vj.  For  convenience,  we  also  define  f  =  0u  and 

f7  =  ^u-  n,  (10) 


as  the  upwinded,  normal  flux  along  the  boundary  of  £. 
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Quadrature-free  approximation:  DG  methods  typically  require  more  degrees  of  free¬ 
dom  than  their  standard  C°  Galerkin  counterparts  and  also  require  the  evaluation  of  both  volume 
and  surface  integrals  over  each  element.  These  factors  may  lead  to  increased  operation  counts 
and  storage  requirements  for  DG  methods,  particularly  when  numerical  quadrature  is  used  to  eval¬ 
uate  element  integrals.  On  the  other  hand,  a  quadrature-free  implementation,  in  which  element 
and  boundary  integrals  are  evaluated  analytically,  can  reduce  operation  counts  and  storage  require¬ 
ments  for  a  DG  discretization  of  Equation  8  significantly  (Atkins  and  Shu  1996).  Although  the 
standard  approach  for  DG  methods  is  to  use  numerical  integration  (Cockburn  and  Shu  1998),  a 
quadrature-free  implementation  for  Equation  8  is  straightforward.  In  essence,  only  two  basic  mod¬ 
ifications  are  necessary  to  accommodate  the  quadrature-free  approach.  The  first  is  a  projection  of  f 
in  [Vh(S)]d,  and  the  other  is  a  representation  of  numerical  fluxes  on  element  boundaries  (Marchan- 
dise  et  al.  2006). 

Before  describing  this  quadrature-free  RKDG  approach,  we  introduce  some  additional  notation. 
An  element  boundary  is  written  as  e,  and  the  global  set  of  element  boundaries  in  M.h  is  {  e,  }.  That 
is,  {e*}  is  the  set  of  vertices  in  A4h  in  one  spatial  dimension,  the  set  of  edges  in  A4h  for  two- 
dimensional  problems,  and  faces  in  three  dimensions.  A  tilde  (~)  is  used  to  distinguish  quantities 
associated  with  element  boundaries. 

In  the  following,  we  restrict  ourselves  to  conforming  meshes  and  introduce  a  local  trial  space  1 4(e) 
on  each  Md_1  element  boundary,  e,  contained  in  A4h  (Atkins  and  Shu  1996).  For  lack  of  better 
notation,  we  label  the  basis  for  14(e)  as  {i\4},  for  i  =  1, . . . ,  hp,  and  the  Md_1  reference  simplex 
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as  e  (e  is  a  point  for  d  =  1).  We  define  Ni  as  before 

Ni  =  fii  o  G  1 

where  G  is  a  consistently  chosen  map  from  e  to  e  and  {iV^},  i  —  1, . . . ,  hp  is  a  basis  for  Pk(e).  A 
natural  fit  is  to  use  a  Lagrangian  basis  for  Pk(e)  with  associated  nodes  p j  E  Md_1. 

For  every  element,  we  require  a  mapping  from  the  elemental  space  Vh(S)  to  14(e)  for  each  e  E  ()£. 
Accordingly,  we  introduce  the  trace  mapping  Tf  :  14(f)  — >  14 (e/) 

Tip 

Tfm  =  J2t^n,  (id 

i= 1 

Ti%  =  Nj(Pi),  Pi  =  G(Pt)  (12) 

In  other  words,  for  each  boundary  face,  ei,  of  an  element  f ,  we  have  a  local  matrix  Tf  e  Mr'?'xn?' 
that  takes  members  of  the  element  trial/test  space  to  a  d—1  dimensional  space  defined  on  the 
boundary  face  that  is  independent  of  the  element’s  local  coordinate  system.  Similarly,  the  transpose 
of  T f  can  be  used  to  map  from  Vh{ei)  to  14(f). 

To  define  ih  e  [14  (f  )]d,  we  use  either  a  coordinate- wise  L2  projection 


fkNi  dx=  4>h,ukNi  dx,  i  =  1, . . . ,  rip,  k  =  l,...,d 


(13) 


or  simple  nodal  interpolation 

fh  J  =  0J«fc(P j),  3  =  1,  •  •  • ,  k  =  1, . . . ,  d  (14) 

A  unique  upwinded  numerical  flux  fup  e  14(e)  is  also  necessary  for  each  e  6  A4h.  We  set 


rip 

fr  =  Y,f:rJN, 

3= 1 


(15) 


with  =  oupdn  ■  ne(p?).  The  unit  outer  noimal,  ne,  is  arbitrarily  chosen  to  point  from  “left” 
to  “right,”  so  that  it  will  be  ±  the  unit  outer  normal  for  the  elements  neighboring  e.  The  value  of 
(j>yk  is  defined  using  the  left  and  right  traces 


fiupj 


= 0r(p  j) 


S(Pj)  u-ne>0 
(j>h(Pj)  otherwise 


(16) 


The  left  and  right  traces  are  defined  simply  by  mapping  the  local  element  degrees  of  freedom, 
<p£  G  Wp,  from  the  neighboring  elements  (f L  and  SR),  to  the  corresponding  degrees  of  freedom 
for  14(e) 

0L  =  T \$L 

=  T  i4>r  (17) 
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~  L  /  R 

where  <fi  '  G  M'",p .  The  subscripts  /  and  l1  are  local  numberings  for  e  for  the  left  and  right 
neighboring  elements,  respectively. 


Inserting  the  new  notation  into  Equation  8,  we  have 


NiNj  da; 


n  d~\~  1  p 

/  f„-  ViV.da;-^  / 

JS  l=1  Jet 


Niffi  ds 


(18) 


for  i  =  1 , ,np  and  f^pt  =  /“pne(  ■  nf.  The  three  spatial  integrals  in  Equation  18  can  now 
be  calculated  analytically  over  each  element,  and  the  corresponding  semi-discrete  system  can  be 
integrated  in  time  using  a  method  of  lines  approach  and  a  class  of  Runge-Kutta  time  discretiza¬ 
tions  (Marchandise  et  al.  2006).  Calculation  of  the  spatial  integrals  for  two-dimensional  triangular 
meshes  is  detailed  in  Appendix  A. 


Time  integration:  To  integrate  Equation  18,  we  use  a  class  of  explicit,  strong  stability  pre¬ 
serving  (SSP)  Runge-Kutta  schemes.  The  SSP  property  is  based  on  the  assumption  that,  for  a 
given  problem,  a  forward  Euler  time  discretization  is  stable  under  a  given  norm  (and  suitable  time 
step  constraint).  An  SSP  method  is  then  one  that  maintains  this  stability  with  possibly  different 
approximation  order  and  time  step  constraint  (Gottlieb  et  al.  2001).  For  our  purposes,  we  can  con¬ 
sider  linear  ordinary  differential  equation  (ODE)  systems,  since  the  mass  matrix  in  Equation  18  is 
invertible  and  limiting  is  not  used  (Cockbum  and  Shu  2001,  Marchandise  et  al.  2006).  That  is,  we 
have  an  ODE  system 


dy 

dt 


Ly 


(19) 


An  s-stage  family  of  Runge-Kutta  methods  that  are  SSP  for  Equation  19  can  be  written 


ym  =  ym_1  _|_  Afn+1Lym_1,  for  m  =  1, . . . ,  s  —  1 

s— 2 

ys  =  ^2  as,kYk  +  aais- 1  (ys_1  +  A tn+1  Ly5-1) 

k= 0 

yn+1  =  y*  (20) 


for  suitably  chosen  coefficients  aStk,  k  —  0, . . . ,  s  —  1  (Gottlieb  et  al.  2001).  A  DG  discretization 
of  order  p  should  be  stable  for  Equation  20  with  the  Courant-Friedrichs-Lewy  (CFL)  condition 


At  < 


h 

c(2p  +  1) 


(21) 


where  h  is  the  element  size  and  c  >  ||u||  for  all  £  G  A4h  (Cockbum  and  Shu  2001,  Marchandise 
et  al.  2006).  For  the  numerical  experiments  below,  we  use  the  k  +  1th  version  of  Equation  20  along 
with  Equation  21  when  the  spatial  approximation  is  based  on  Pk  elements.  The  corresponding 
coefficients  can  be  found  in  Table  3.1  of  Gottlieb  et  al.  (2001). 


NUMERICAL  EXPERIMENTS:  To  evaluate  our  quadrature-free  RKDG  approach,  we 
consider  several  classical  test  problems  for  propagating  interfaces  with  a  specified  velocity  (Rider 
and  Kothe  1995,  Sussman  and  Fatemi  1999,  Pilliod  and  Puckett  2004,  Olsson  and  Kreiss  2005, 
Marchandise  et  al.  2006). 
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Test  problems:  The  physical  domain  for  the  first  test  problem  (PA)  is  0  =  [0, 1]  x  [0,1], 
and  the  velocity  field  is  u  =  [27 r(y  —  1/2),  27t(1/2  —  x)].  The  exact  solution  is  a  cone  of  radius 
r0  =  1/8.  A  parameterization  for  the  exact  solution,  oex,  is 


(pex(x,y,t) 


(1  +  cos(7rx/r0))  (1  +  cos(ny/r0))  /4,  ||r||  <  r0 
0  otherwise 


(22) 


where  x  =  x  —  xc,  y  =  y  —  yc,  and 

xc  =  sin(27rf)/4  +  1/2,  yc  =  cos(27rf)/4  +  1/2 


The  second  problem  (PB)  is  set  on  Q  =  [0, 1]  x  [0, 1]  with  a  velocity  field  given  by 

ux  =  cos(7rf/8)  sin(27T7/)  sin2(7rx) 

uy  =  —  cos(7rf/8)  sin(27rx)  sin2(7T7/)  (23) 

The  initial  condition  is  a  disk  of  radius  0.15  centered  at  (0.5,  0.75),  with  an  initial  signed  distance 
function  d0  =  (x  —  0.5)2  +  (y  —  0.75)2  —  0.152.  The  solution  should  return  to  the  initial  condition 
at  T  =  8. 

The  final  example  (PC)  is  on  =  [0, 100]  x  [0, 100]  with  a  velocity  field 

u  =  (7t(50  —  y)/ 314, 7r(x  —  50)/314)T 

The  initial  condition  is  Zalesak’s  slotted  disk  with  a  radius  of  15,  width  of  5,  and  slot  length  of  15 
(Sussman  and  Fatemi  1999,  Marchandise  et  al.  2006). 

Illustrative  results:  For  each  problem,  a  regular  triangulation  of  Vt  was  formed  with  Nx 
triangles  along  the  x  axis,  and  Ny  triangles  along  the  y  axis.  Simulations  were  performed  with 
varying  orders  of  approximation  from  k  =  1, . . . ,  4.  A  target  Courant  number,  Cr,  was  chosen 
close  to  the  maximum  allowed  for  each  k  as  given  by  Equation  21.  A  dual  processor  Mac  G5  (2 
GHz)  with  2  GB  RAM  was  used  for  the  computations.  The  methods  were  implemented  in  C+  + 
and  compiled  with  gcc  version  3.3  and  -O  optimization. 

Table  1  summarizes  the  computations  performed  for  PA,  while  Table  2  contains  the  corresponding 
L 1  errors,  CPU  times,  and  mass  errors.  The  mass  error  here  is  simply  defined  to  be  the  normalized 
difference  between  the  mass  in  the  numerical  and  analytical  solutions.  In  Table  2,  Nj  is  the  total 
number  of  degrees  of  freedom  on  each  mesh.  Figure  3  shows  the  initial  condition  and  solution  at 
T  =  0.5  for  a  P 2  approximation  on  a  grid  with  64  x  64  triangles. 

Table  3  summarizes  the  computations  performed  for  PB,  and  Table  4  contains  the  corresponding  L1 
errors,  total  mass  errors,  and  CPU  times.  Figure  4  shows  the  zero  contour  of  the  numerical  solution 
for  a  P 2  approximation  on  a  64  x  64  regular  grid.  The  exact  solution  for  PB  is  identical  to  the  initial 
condition.  The  right  plot  illustrates  the  significant  deformation  that  the  solution  experiences  before 
returning  to  the  initial  state  at  T  =  8.  Note  the  difference  in  scales  for  the  two  plots. 

The  simulations  performed  for  PC  are  summarized  in  Table  5,  while  Figure  5  illustrates  the  rel¬ 
ative  accuracy  for  first-  through  fourth-order  approximations  on  the  same  128  x  128  grid.  The 
difference  in  accuracy  is  most  evident  around  corners,  where  the  lower  order  methods  are  smeared 
significantly  (Sussman  and  Fatemi  1999,  Marchandise  et  al.  2006). 
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Table  1 

Simulation  summary  PA 

Run  Problem 

pk 

Nx 

Ny 

Nd 

Cr 

PA.l 

PA 

1 

128 

128 

98304 

0.3 

PA. 2 

PA 

2 

64 

64 

49152 

0.18 

PA. 3 

PA 

3 

32 

32 

20480 

0.128 

PA.4 

PA 

4 

16 

16 

7680 

0.1 

Table  2.  Simulation  results  PA 

Run 

Li  Error 

Mass  Error 

CPU  [s] 

PA.l 

1.3097  x  10"4 

0.000 

3.214  x  102 

PA. 2 

1.1110  x  10~4 

0.000 

2.048  x  102 

PA. 3 

1.5038  x  10”4 

5.600  x  10~6 

8.237  x  101 

PA.4 

3.3699  x  10“4 

2.230  x  10"5 

2.805  x  101 

Table  3.  Simulation  summary  PB 

Run 

Problem 

pk 

Nx 

Ny 

Cr 

PB.l 

PB 

1 

128 

128 

0.3 

PB.2 

PB 

2 

64 

64 

0.18 

PB.3 

PB 

3 

64 

64 

0.128 

PB.4 

PB 

4 

32 

32 

0.1 

DISCUSSION:  The  numerical  results  above  are  preliminary.  They  indicate  that  a  quadrature- 
free  RKDG  approach  can  accurately  resolve  propagating  interfaces  for  simple  rotating  velocity 
fields  and  more  complex  two-dimensional  flows  as  in  PB  (Rider  and  Kothe  1995).  The  results 
show  that  higher  order  approximations  are  more  accurate  for  the  same  spatial  resolution.  In  some 
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[left].  0  at  T  =  1.8  [right] 


Table  4.  Simulation  results  PB 

Run 

Li  Error 

Mass  Error 

CPU  [s] 

PB.l 

2.6219  x  10~3 

0.000 

1.416  x  103 

PB.2 

1.6250  x  10"3 

0.000 

8.199  x  102 

PB.3 

1.0616  x  10"3 

0.000 

2.566  x  103 

PB.4 

1.6091  x  10"3 

1.000  x  10~6 

8.113  x  102 

Table  5.  Simulation  summary  PC 

Run 

Problem 

pk 

Nx 

Ny 

Cr 

PC.l 

PC 

1 

128 

128 

0.3 

PC. 2 

PC 

2 

128 

128 

0.18 

PC. 3 

PC 

3 

128 

128 

0.128 

PC. 4 

PC 

4 

128 

128 

0.1 

Figure  5.  PC,  T  =  628.  Zero  contour,  P 1  (blue),  P 2  (red),  P3  (green),  and  P4  (black) 
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cases,  similar  accuracy  can  be  achieved  for  much  coarser  meshes  with  higher  order  approximations 
and  an  overall  reduction  in  computational  effort.  However,  this  work  does  not  constitute  a  rigorous 
evaluation  of  computational  efficiency.  Moreover,  the  methods’  performance  needs  to  be  evaluated 
for  problems  in  which  there  is  nontrivial  coupling  between  the  interface  location  and  fluid  flow. 

ADDITIONAL  INFORMATION:  This  CHETN  is  a  product  of  the  High  Fidelity  Vessel  Ef¬ 
fects  Work  Unit  of  the  Navigation  Systems  Research  Program  being  conducted  at  the  U.S.  Army 
Engineer  Research  and  Development  Center,  Coastal  and  Hydraulics  Laboratory.  Questions  about 
this  technical  note  can  be  addressed  to  Dr.  Christopher  Kees  (Voice:  601-634-3110, 
email:  Christopher  .  e  .  kees@erdc  .  usace  .  army  .mil).  For  information  about  the  Nav¬ 
igation  Systems  Research  Program,  please  contact  the  Navigation  Systems  Program  Manager, 
Dr.  John  Hite  (Voice:  601-634-2402,  email:  John  .  E  .  Hite@erdc  .  usace  .  army  . mil).  This 
technical  note  should  be  cited  as  follows: 

Farthing,  M.,  and  C.  Kees.  2008.  Implementation  of  discontinuous  Galerkin  methods 
for  the  level  set  equation  on  unstructured  meshes.  Coastal  and  Hydraulics  Engineering 
Technical  Note  ERDC/CHL  CHETN  XIII-2.  Vicksburg,  MS:  U.S.  Army  Engineer 
Research  and  Development  Center. 
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APPENDIX  A.  ELEMENT  INTEGRAL  CALCULATIONS:  To  illustrate  the  cal¬ 
culations  necessary  to  evaluate  Equation  18  analytically,  we  consider  a  two-dimensional  example. 
The  element  matrix  on  the  left-hand  side  of  Equation  18  is  the  standard  mass  matrix 


Mij  =  j  NiNj  dx  =  J  NiNj\  det  J\  dxdy  =  |  det ./ 1 M%3  (24) 

where  J  is  the  Jacobian  of  F.  The  value  of  J  is  constant  on  each  element,  since  F  is  affine,  and 
M  =  (Mij)  €  MnPxnp  is  only  a  function  of  the  reference  element  and  the  local  shape  functions 
chosen. 


The  first  integral  on  the  right-hand  side  of  Equation  18  is  also  straightforward  to  calculate.  Recall- 


vNi  and  setting 

=  [fh  fHY » we  have 

j  f  ■  ViVi  dx  = 

Up 

|detJ|^/oz3i  + 
3= 1 

Up 

|detJ|^s"D| 

3= 1 

(25) 

g*’j  = 

fX'j  J 11  + 

(26) 

gy,j  = 

+  fy,3j-i 

(27) 

Di  = 
v 

i  dxNldxdV 

(28) 

D*.  = 

1-3 

l 

(29) 

Using  the  local  trace  mapping  (transposed)  to  write 

Up 

N<  =  J2Ti£iA  (30) 

3= 1 


the  boundary  integral  in  Equation  1 8  can  be  written  in  terms  of  the  element  boundary  flux  degrees 
of  freedom,  for  j  =  1 , ,hp  and  local  element  matrices 


«f  ki 


E  ey, 


up, 3 
ei 


3= 1 


V  Mivlfj 


jk 


k= 1 


(31) 

(32) 


Ml,ij  — 


NiNj  ds 


'  et 


(33) 


Here,  af  =  ng  ■  n,:(  accounts  for  the  current  element’s  outer  normal  orientation  relative  to  the 
unique  global  orientation  assigned  to  the  boundary  face. 
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For  each  £  e  Aih,  we  can  define  element  degree  of  freedom  vectors  cj)£,  gf,  and  g|  in  Mnp  and 
for  each  element  boundary,  e*,  a  vector  of  numerical  flux  degrees  of  freedom,  f“p  e  Mnp.  We  then 
have  the  local  elemental  equation 


r)fh  d-\-l  |  | 

=  D'gf  +  D»g|  -  af  ^  (34) 

1=1  I  I 

APPENDIX  B.  P 2  LOCAL  MATRICES  FOR  TRIANGLES:  Finally,  we  repro¬ 
duce  the  local  element  matrices  for  a  triangular  mesh  with  a  P 2  local  approximation  space.  The 
barycentric  coordinates  on  £  are 

A0  =  1-x-y 
Xi  =  x 

A  2  =  y  (35) 

The  P2  shape  functions  are  then 

=  A0(2A0  —  1) 

N,  =  A1(2A1-1) 

-^2  —  A2  (2A2  —  1) 

-^3  =  4A0Ai 

N4.  =  4A1A2 

N5  =  4A0A2  (36) 

The  nodal  locations  associated  with  {iV*},  %  —  0, . . . ,  5  are  illustrated  in  Figure  6.  Note  that,  in  the 
matrix  definitions  below,  the  edges  are  numbered  counterclockwise  around  £  starting  with  the  x 
axis. 
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